# In this script we perform:

#  1. full asset-liability matching using 15 bonds, with maturities spanning over
#     each of the 15 time periods - the ones used in script
#     "S3.Estimating the market yield curve by replication.R".

#  2. cash flow immunization using 3 bonds, defined below.

library(tidyverse)
library(readxl)

############### Example of Asset-Liability Matching #############

# Load the market yield curve and the matrix of payments of bonds, taken from
# the same script.

load("Examples/market.yield.curve.Rdata")
market.yield.curve

payments.tbl <- read_xlsx(
  "Examples/interestRateData2025.xlsx", sheet = "paymentMatrixC"
)

# Consider matrix C
C 




# Assume a liability cash flow of 1M per time period.
data <- market.yield.curve 




# Let us find the weights of the (full asset-liability) matching portfolio:
# (w_1,...,w_n) = (L_1,...,L_n) * inverse(C)
w.match




# Verification


#################### Example of Cash Flow Immunization ################

# Most times, cash flow matching is not practical.
# We will perform cash flow immunization using 3 bonds.
maturities <- c(2, 8, 11)
coupons <- c(0.04, 0.05, 0.05)
face.values <- c(100, 100, 100)

# Let us add the bond payments to our data
data <- data %>%
  # Include the coupon payments
  mutate(
    coupons.1 = ,
    coupons.2 = ,
    coupons.3 = 
  ) %>% 
  # Include the principal payments
  mutate(
    principal.1 = ,
    principal.2 = ,
    principal.3 = 
  ) %>%
  # Compute total payments
  mutate(
    payments.1 = ,
    payments.2 = ,
    payments.3 = 
  ) %>% 
  # Select the columns that matter
  select()


# Compute the duration and convexity of the liabilities cash flow (in the same 
# way as in script "S4.Estimating Present Value Sensitivity to Yield Changes.R")
sens.liab <- data %>% 
  mutate(PV = ) %>% 
  mutate(
    d_aux = ,
    c_aux = 
  ) %>% 
  summarize(
    PV = ,
    d_aux = ,
    c_aux = 
  ) %>% 
  mutate(duration = , convexity = ) %>%
  # Select the columns that matter
  select()
sens.liab

# Now we compute the duration and convexity of each of the 3 bonds
# We will encapsulate our code in a function in order to avoid repetition when 
# computing duration and convexity of our three bonds. Please note that when
# we invoke a column from a data frame inside a function using tidyverse, we
# have to "embrace" it, meaning that we surround it with {{...}}.

f.sensitivity <- function(table, variable){
  
  table 
  
}

sens.bond.1 <- f.sensitivity(table = , variable = )
sens.bond.2 <- f.sensitivity()
sens.bond.3 <- f.sensitivity()
sens.bond.1
sens.bond.2
sens.bond.3

# In order to find the correct amounts of each bond to include in our portfolio W,
# we will construct the system of equations prescribed for second order immunization:
# PV_L = PV_A                             (equate present values)
# PV_L * Duration_L = PV_A * Duration_A   (equate dollar durations)
# PV_L * Convexity_L = PV_A * Convexity_A (equate dollar convexities)

# Function for computing Dollar Duration and Dollar Convexity
f.dollar <- function(sens.tbl){
  
  sens.tbl %>%
    mutate(DD =, DC = ) %>% 
    select()
  
}

PV.sens.liab <- f.dollar(sens.liab)
PV.sens.liab
PV.sens.bond.1 <- f.dollar(sens.bond.1)
PV.sens.bond.1
PV.sens.bond.2 <- f.dollar(sens.bond.2)
PV.sens.bond.2
PV.sens.bond.3 <- f.dollar(sens.bond.3)
PV.sens.bond.3

# Our system of equations is of the form M * w = v, where:
#  - M is a matrix corresponding to the PV and sensitivity of the assets (known)
#  - w is the vector of portfolio weights (unknown)
#  - v is a vector corresponding to the PV and sensitivity of the liabilities (known)

# t() is the 'transpose' function in base R. Transpose the three vectors (PV, DD, DC)
# we have obtained above and bind them in a matrix (M)
t(PV.sens.bond.1)
t(PV.sens.bond.2)
t(PV.sens.bond.3)
M <- cbind(t(PV.sens.bond.1), t(PV.sens.bond.2), t(PV.sens.bond.3))
M

# What is the vector of our objectives for PV, DD and DC?

v

# Solving the system of equations (w = inverse(M) * v):

w

# These weights are our solution. We have determined an immunizing portfolio to
# our liabilities. Let us check the asset cash flow, as well as its sensitivity
# measures.
data
data <- data %>% 
  mutate(
    # Cash flows for each bond
    cf.1 = ,
    cf.2 = ,
    cf.3 = ,
    # Assets cash flow
    assets.cf = cf.1 + cf.2 + cf.3
  )

# Same D and C as liabilities? Please check!
data 

# Same DD and DC as liabilities? Please check!
data %>% f.sensitivity(assets.cf) %>% f.dollar() 

# Calculate the present value of the surplus
data <- data %>%
  mutate(
    surplus =  assets.cf - liab.cf,
    PV.liab = P_t * liab.cf,
    PV.assets = P_t * assets.cf,
    PV.surplus = P_t * surplus
    )

data %>% summarize(sum(PV.surplus))

# Save the immunization data for later
immuniz.data <- data %>% 
  select(t, liab.cf, assets.cf, surplus, starts_with("PV"))
save(immuniz.data, file = "Examples/immuniz.data.Rdata")

totals <- data %>%
  summarize(
    across(c(liab.cf, assets.cf, surplus, PV.liab, PV.assets, PV.surplus), sum)
    )
# The PV of our surplus is approximately zero
totals

# Finally, we will see how our portfolio holds if we add a few perturbations
# to our yield curve. We include a perturbation 'a' for the first three years
# and a perturbation 'b' for the remaining years.
a <- 0.005
b <- -0.005
yield.perturbation <- c(rep(a,3),rep(b,12))
yield.perturbation

# Homework: compute the present value of the liabilities, assets and surplus in
# this new case. It is easy!

perturbed.data <- data




# PV of assets, PV of liabilities and PV of surplus
perturbed.data






# Conclusion: